library(spatstat)
library(rgdal)
library(maptools)
data(redwoodfull)
d <- redwoodfull

##### read from shapefile --not working!
hs <- readOGR("/home/perry/data/hotsprings/hot_springs.shp","hot_springs")
hs.ca <- hs[hs$SC == "CA",]
crds <- coordinates(hs.ca)

#states <- readOGR("/home/perry/data/states/statesp020.shp","statesp020")
ca <- states[states$STATE == "California",]
bbox <- summary(ca)$bbox
#d <- as(hs.ca,"ppp")
xrange <- bbox[3]-bbox[1]
yrange <- bbox[4]-bbox[2]
xx <- xrange*0.025
yy <- yrange*0.025

d <- ppp(crds[,1], crds[,2], c(bbox[1]-xx,bbox[3]+xx), c(bbox[2]-xx,bbox[4]+xx))

pdf("~/Desktop/spatstat_test.pdf", width=17,height=11)

#######
a <- plot(d, pch=16, cex=0.5, main="Quadrat Map")
q <- quadratcount(d)
plot(q, add=T)

plot(ca, add=T)


#######
# If the point pattern has "marks" i.e. a categorical attribute
#legend(-0.25, 0.5, names(a), pch = a)
#plot(split(d))

#### Find intensity of point pattern
lamb <- summary(d)$intensity

#######
oldpar <- par(mfrow=c(1,2))
k <- Kest(d)
g <- Gest(d)
plot(k)
plot(g)

#######
par <- oldpar
dens <- density(d)
plot(ca,col="green",main="Intensity of Density")
#plot(dens, main="Density")
points(d,pch=16,cex=0.25) 
contour(dens, add=T); 

persp(dens, main="Density")

#sapply(seq(0.03,0.11,0.01), 
#  function (x) { 
#    plot(density(d,x), main=paste("Density with sd=",x))
#    points(d,pch=16,cex=0.25) 
#    contour(density(d,x), add=T); 
#  } 
#)

###### 3D plot
#par <- oldpar
#oldpar <- par(mfrow=c(3,3))
#sapply(seq(0.03,0.11,0.01), 
#  function (x) { 
#    persp(density(d,x), main=paste("Density with sd=",x))
#  } 
#)

dev.off()
